%%Multiple images
function [Z,R] = multiple_images(path,C,S,WE)
images = dir([path '/*.jpg']);
num_img = length(images);
for img=1:num_img;
    I(img,:,:,:) = cast(imread([path '\' images(img).name]), 'double')./255;
end

%% get weight maps for all images
for img=1:num_img;
    F(img,:,:,:) = combined_weights(squeeze(I(img,:,:,:)),C,S,WE) + 1e-12;
end

%% normalize weight maps
[num_img,height,width,rgb] = size(I);

for img=1:num_img;
    N(:,:,img) = zeros(height,width);
    Den = sum(F);
    N(:,:,img) = squeeze(F(img,:,:)./Den);
end

%%
Z = zeros(height, width, rgb);
for img=1:num_img;
    for i = 1:height
        for j = 1:width
            for k = 1:rgb
                Z(i,j,k) = Z(i,j,k) + I(img,i,j,k).*N(i,j,img);
            end
        end
    end
end

figure, imshow(Z), title('Naive Blending', 'FontWeight','bold')

%Initalize empty pyramid with height/width of images to get maximum number
%of praymid levels
pyramid = gauss_pyramid(zeros(height,width,3),height,width);
numL = length(pyramid);

%Build the pyramids for each image
for i = 1:num_img
	pyramidW = gauss_pyramid(N(:,:,i),height,width); %Weight map pyramid
	pyramidI = laplace_pyramid(squeeze(I(i,:,:,:)),height,width); %Image pyramid
    
    %Blend images using weight maps at each pyramid level
    for l = 1:numL
        for j = 1:3
            weight_map(:,:,j) = pyramidW{l}; %Use same weight map for all colour channels
        end
        pyramid{l} = pyramid{l} + weight_map.*pyramidI{l}; %Blend images using weight map
        clear weight_map; %weight_map changes size from level to level, need to reset variable
    end
end

%Rebuild image from blended pyramid
R = recon_laplace_pyramid(pyramid);
figure, imshow(R), title('Laplacian Pyramid / Multiresolution Blending', 'FontWeight','bold');